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We devise a soft, solvent-free, coarse-grained model for lipid bilayer membranes. The non-bonded 
interactions take the form of a weighted-density functional which allows us to describe the ther- 
modynamics of self-assembly and packing effects of the coarse-grained beads in terms of a density 
expansion of the equation of state and the weighting functions that regularize the microscopic bead 
densities, respectively. Identifying the length and energy scales via the bilayer thickness and the 
thermal energy scale, fcsT, the model qualitatively reproduces key characteristics (e.g., bending 
rigidity, area per lipid molecules, and compressibility) of lipid membranes. We employ this model 
to study the main phase transition between the liquid and the gel phase of the bilayer membrane. 
We accurately locate the phase coexistence using free energy calculations and also obtain estimates 
for the bare and the thermodynamic line tension. 
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I. INTRODUCTION 

Lipid bilayers are one of nature's most ingenious 
inventions*^ They serve as a compartment to all cells, 
which form the building blocks of life, and they mediate 
the transport of molecules from the inside to the outside 
of cells. Many important properties of bilayer membranes 
involve collective phenomena, where a large number of in- 
teracting lipid molecules participate. Examples include 
the self-assembly of amphiphilic molecules into bilayer 
membranes, phase transitions between different phases or 
changes of the membrane topology, e.g., pore formation 
or fusioni^ Computer simulations contribute to the un- 
derstanding how these collective phenomena depend on 
the properties of the individual, constituent molecules. 

Often, collective phenomena involve mesoscopic time- 
and length scales - microseconds and nm - which are dif- 
ficult to observe directly in experiments and which are at 
present beyond the scales that can be addressed by mod- 
els with atomistic resolution. Therefore, several compu- 
tational models have been developed where local, atom- 
istic degrees of freedom have been integrated out 4^ The 
reduced number of degrees of freedom of coarse-grained 
models and the softer interactions between the effective 
interaction centers opens up the opportunity to computa- 
tionally address the mesoscopic scales involved in collec- 
tive phenomena in complex biological matter. In order to 
design a coarse-grained model, first, one has to decide, 
which are the relevant degrees of freedom for the phe- 
nomena under study and which are to be integrated out. 
Second, one has to construct the effective interactions 
between the remaining degrees of freedom. This con- 
struction is either performed systematically by explicitly 
tracing out the microscopic degrees of the freedom, or 
one invokes the concept of universality and uses a min- 
imal set of interactions that is comprised only of those 
interactions, which are necessary to bring about the phe- 
nomena under study. The strength of those relevant in- 
teractions can be parameterized by comparing properties 



of the model to experimental data. 

In this study we rely on the concept of universality. 
Instead of trying to reproduce chemical details of a spe- 
cific lipid molecule, like it is done in atomistic or system- 
atically coarse-grained simulations, we present a coarse- 
grained, solvent-free model for amphiphilic bilayers. Our 
model has some similarities with models used in self- 
consistent field calculations and it interpolates smoothly 
between lipid bilayers and polymeric membranes. Within 
the mean-field approximation of our model, there is a 
clear separation between the thermodynamic properties 
and the local fluid structure (i.e., packing effects) of the 
hydrophobic core of the bilayer membrane. We inves- 
tigate the main phase transition between a liquid and 
a gel state of a self-assembled, one-component bilayer 
membrane and the line tension between domains. The 
main phase transition has been characterized in exper- 
iments for many lipids; 7 -^ and it has also been consid- 
ered in many coarse-grained models of lipid bilayers^r— 
The relation between the microscopic properties of the 
lipid molecules (e.g., the stiffness of the hydrocarbon 
tails and the fluid-like packing effects) and the macro- 
scopic phase behavior, however, is only incompletely un- 
derstood. Moreover, there are only very few attempts to 
accurately locate the phase boundaries. 17 Problems arise 
from the hysteresis effects and metastability at the first- 
order transition, which seriously hamper the accurate de- 
termination of the location of the main phase transition 
by computer simulation. Additionally, the line tension 
between the fluid and gel phases has not been measured 
in coarse-grained models, which retain the notion of lipid 
molecules. This free energy of the domain boundaries be- 
tween laterally coexisting phases has attracted abiding 
experimental interest ^Sr— 

In order to describe the main phase transition, our 
coarse-grained model has to incorporate both (i) the 
minimal interactions that bring about the self-assembly 
into a bilayer membrane and, additionally, (ii) further 
details of the local inter- and intramolecular structure 
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that give rise to the transition from a liquid to a gel 
phase. Thus, the following relevant properties are re- 
tained in our coarse-grained representation: (i) Each lipid 
comprises two different constituents, a hydrophobic tail 
and a slightly smaller hydrophilic head, which repel each 
other. These interactions drive the self-assembly into bi- 
layer membranes. Since the lipid bilayer is typically sur- 
rounded by a solvent, the hydrophilic heads turn towards 
the solvent and the hydrophobic tails lump together in 
the bilayer's center, (ii) The hydrocarbon tails of the 
lipids are characterized by a finite length, a limited con- 
formational flexibility and a finite excluded volume diam- 
eter. These interactions give rise to a crystalline packing 
of the molecules in the gel phase. Incorporating both 
aspects, our coarse-grained model bridges between min- 
imal coarse-grained representations, which only capture 
the universal aspects of self-assembly, and systematically 
coarse-grained models that have been explicitly derived 
from an atomistic model. 

Our manuscript is arranged as follows: In Sec. |TT] 
we describe our solvent-free, coarse-grained model with 
soft interactions and provide details of the Multibody 
Dissipative Particle Dynamics (MDPD) simulations. A 
technical description of the symplectic integration algo- 
rithm for simulating in an ensemble with constant ten- 
sion is deferred to Appendix [A] The subsequent section, 
Sec. IIII1 demonstrates the self-assembly of lipids into bi- 
layer membranes. Several static and dynamic proper- 
ties of our model, such as the bending rigidity, the bi- 
layer density profile, and the molecular diffusion coeffi- 
cient are measured. In Sec. IIVI the main phase transi- 
tion is studied. We use Umbrella Sampling (US) 2 ^— 
to restrain the fluctuations of an orientational order pa- 
rameter. Changing the order parameter, we reversibly 
transform the fluid into a gel phase and obtain the con- 
comitant free-energy profile by the Weighted Histogram 
Analysis Method (WHAM)^.^^ The bilayer configu- 
rations along the reversible path are discussed. Sec.lVlde- 
scribes the measurement of the thermodynamic line ten- 
sion between gel and fluid domains, extracted from the 
free-energy profile, and the bare line tension, computed 
from the fluctuation spectra of the domain boundaries. 
The relation between these two properties is discussed in 
Appendix |Bj The paper concludes with a summary and 
an outlook in Sec. I VII 



II. MODEL AND TECHNIQUE 
A. Model 

We consider a coarse-grained model for the simula- 
tion of lipid bilayer membranes. Our system contains n 
lipid molecules that are represented by linear bead-spring 
chains comprising N — 16 effective interaction cen- 
ters, which are either hydrophobic ("A") or hydrophilic 
("B"). The ratio of A-beads, Na, in a lipid is defined 
by the asymmetry parameter /, so that Na = fN and 



Nb = (1 — f)N, respectively. The beads are connected 
by harmonic springs with spring constant k s . Addition- 
ally, we apply a bond-angle potential between every three 
successive beads i — i + 1 with constant k\, to stiffen 
the lipids. Thus, the intramolecular, bonded interactions 
of a single lipid are given by 

14 N ~ x b N ~ 1 
■ Ji '-- - V- Ik [rm _ ri f + J2 h [1 - cosOi] , (I) 



k B T 
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i=l 



i=2 



where 9i is the angle between the vectors (r^ — r^_i) 
and (rj_)_i— rj). The thermal energy, k B T, serves 
as the unit of energy in our model. Although lipid 
molecules are characterized by several length scales, we 
use the root-mean-squared end-to-end distance, R eo = 
((rj — r^v) 2 ) 1 / 2 , of lipids that are only subjected to the 
bonded interactions, %b ) 30 ' 31 as the characteristic dimen- 
sion of the bilayer. It can be pictured as the head-to-tail 
length of a single lipid in vacuum. The use of R eo to 
specify the molecular extension is rooted in polymeric 
membranes, where the polymer conformations are char- 
acterized by this single length scaled Its value, in turn, 
depends on the values of N, k s , and kb- The bond stiff- 
ness restricts the conformational fluctuations of the am- 
phiphilic molecules such that the average molecular size 
and its shape fluctuations are controlled by the parame- 
ters of the model. The actual size of a lipid molecule, of 
course, is influenced by the interactions with its neigh- 
bors, e.g., it differs in the liquid and the gel phase. 

Since on large length scales a bilayer membrane can be 
conceived as a thin, two-dimensional sheet embedded in 
a three-dimensional volume, most of the volume is occu- 
pied by solvent. Although the solvent acts as a transport 
medium in a plethora of biological processes and mediates 
the self-assembly, drastically simplifying its representa- 
tion or even integrating out the solvent altogether offers 
a potentially huge reduction in the number of the degrees 
of freedom i 13 ' 14 i 16 ' 32 " — By integrating out the degrees of 
freedom of the solvent, the original interactions of the un- 
derlying model containing the explicit solvent molecules 
are turned into effective interactions. These depend on 
the thermodynamic state, at which the elimination of 
the explicit solvent has been performed. Thus, the non- 
bonded interactions are free energies and care has to be 
exerted when extracting thermodynamic properties^ 

In the following we employ a solvent-free model to 
study thermodynamic equilibrium properties. Thus, hy- 
drodynamic interactions, that are mediated by the sol- 
vent, are irrelevant. The non-bonded interactions are ac- 
counted for by a phenomenological Ansatz for the excess 
free energy. Specifically, we use an expansion up to third 
order for the non-bonded excess free energy in terms of 
the dimensionless, weighted densities of the molecules^S 
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A summation over all Greek indices that occur twice is 
implied and the integration extends over the whole vol- 
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ume of the simulation box. The term in the bracket de- 
notes the excess free energy per particle. The weighted 
densities are related to the explicit particle coordinates 
via a weighted average over a small volume. The details 
of this procedure are discussed below. Here we only note 
that, once the weighted densities are specified in terms 
of the microscopic particle coordinates, the Hamiltonian 
(12) becomes a function of the explicit particle coordi- 
nates and the properties of the coarse-grained model can 
be studied by computer simulation. 

Within the mean-held approximation, the properties 
of the particle-based simulation model coincide with the 
results of a density functional theory (DFT) calculation 
using the excess free energy functional H' nh [pA, Pb}- In 
particular, within the mean-held approximation, thermo- 
dynamic and structural properties decoupled The ther- 
modynamic properties of a spatially homogeneous sys- 
tem, e.g., the equation of state, are dictated by the seven 
expansion coefficients, v a p and w a ^. The local struc- 
ture of the liquid, in turn, is encoded in the definition of 
the weighted densities. 

The advantages of these DFT-based, non-bonded in- 
teractions are twofold: On the one hand, Eq. <j2j) can be 
generalized in a systematic way to accommodate more 
sophisticated equations of state. In the present work, 
we use a third-order expansio n 30 i 39 because this is the 
simplest form capable of describing all six, qualitatively 
different types of phase diagrams that a compressible bi- 
nary system exhibits accord ing to the classification of 
Ivan Konvnenburg and Scottl 4£ i.e., it suffices to capture 
all qualitative features of the interplay between liquid- 
vapor phase separation and demixing of two species. 
Moreover, by virtue of its simplicity, the second- and 
third-order coefficients are straightforwardly related to 
the density and compressibility of a homogeneous liq- 
uid and the incompatibility between hydrophobic and 
hydrophilic entities. This relation imparts a transpar- 
ent physical interpretation onto the coefficients. The 
density expansion also allows for a systematic general- 
ization to systems comprised of more than two different 
speciesj 41 i 42 This situation naturally arises in the study of 
more complex systems. On the other hand, the weighted 
densities encode local structural information. Altering 
the definition of the weighted density, we are able to de- 
scribe lipid bilayer membranes, which exhibit pronounced 
packing effects on the length scale of an effective interac- 
tion center, or polymersomes that are comprised of long, 
flexible, amphiphilic polymers and, typically, do not form 
gel phases. 

We discuss how to choose the expansion coefficients 
and the definition of the weighted densities in turn. 



1. Thermodynamic coefficients of the third-order density 
expansion 

Formally, we consider the system of amphiphiles and 
solvent on the mesoscopic scale of a coarse-grained in- 



teraction center as an incompressible, dense liquid with 
bulk density p a . Knowing the local densities of am- 
phiphiles, one can reconstruct the solvent density by as- 
suming that the total system of solvent and amphiphiles 
is nearly incompressible and integrate out the degrees of 
freedom associated with the solvent This gives rise 
to effective interactions and the incompressibility con- 
straint generates multi-body interactions. The occur- 
rence of multi-body interactions is natural in the course 
of coarse-graining and it would also arise during a system- 
atic coarse-graining procedure where microscopic degrees 
of freedom are explicitly integrated out. 

The coefficients vaa and waaa dictate the properties 
of the hydrophobic species in contact with the solvent. 
In a solvent-free model, the hydrophobic species forms a 
dense liquid that coexists with a vapor phase, which rep- 
resents the solvent. Since the solubility of amphiphiles in 
the solvent is vanishingly small, the (osmotic) pressure 
of the vapor phase, which coexists with the liquid, van- 
ishes, P w 0. Using the mean-field equation of state for 
the pure ^-component 
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we obtain for the molecular density, p CO cx, of the liquid 
with P = 



Pc 



3VAA 



'{waaa 

and for the dimensionless, inverse compressibility 
R. 
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respectively. In both cases we have neglected the con- 
tribution of the first term in the equation of state © 
that corresponds to an ideal gas. These approximate ex- 
pressions provide a simple physical interpretation of the 
expansion coefficients. We will present our results as a 
function of kN and p COC x using the dependencies 



kN 



VAA 



and 
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WAAA 



3kN + 2 
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(7) 



We use Pcocx as control parameter to study the main 
phase transition between a fluid and a gel phase. At 
large p cocx molecules strongly overlap, packing effects are 
small, and the system is in the fluid phase. This behavior 
is typical for polymersomes, where a coarse-grained bead 
is comprised of many atomistic units or for high tem- 
peratures, where the soft, non-bonded interactions are 
weak compared to the thermal energy scale. A decrease 
of Pcocx: in turn, corresponds to an increase of the repul- 
sive, third-order interactions (cf. Eq. (0), which gives 
rise to a transition from the fluid to the gel phase. 
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The coefficient, vab, sets the strength of the interac- 
tions between A and B beads. It is related to the Flory- 
Huggins parameter, ^JV, via 

\N 1 

WAS = h ~ (^Aj4 + vbb) ■ (8) 

Pcocx ^ 

The dimensionless, invariant quantity, ^TV, measures 
the incompatibility between hydrophilic and hydropho- 
bic species, vbb and wbbb are chosen, such that the 
hydrophilic beads are in a good solvent, i.e., their inter- 
actions are purely repulsive, vbb = 0.1 and wbbb = 0. 
The mixed, third-order coefficients, waab and wabb, do 
not influence the qualitative behavior and, for simplicity, 
we set w A aa = waab = wabb- 

Four phenomenological parameters describe the ther- 
modynamics of our soft, solvent- free, coarse-grained 
model: p CO cx, kN,xN and R eo , which parameterize (i) 
the density and (ii) the limited compressibility of the hy- 
drophobic interior, (iii) the incompatibility between hy- 
drophilic and hydrophobic beads, and (iv) the spatial ex- 
tension of a lipid molecule. All these parameters are di- 
rectly related to experimentally accessible quantities and 
our model can be related to a specific system by match- 
ing these four parameters of our coarse-grained model to 
experimental data. 

For instance, we estimate the order of magnitude of 
kN from the bulk properties of an alkane liquid. Us- 
ing the isothermal compressibility under standard condi- 
tions kt = 0.955 GPa' 1 for n-Dodecane^ its bulk mass 
density p m = 748.8 kg/m 3 , and its molar mass m = 
170.34 g/mol, we obtain nN — m/ (nrPrnkBT) « 98. 

2. Weighted densities 

For lipid bilayer membranes we seek for weighted den- 
sities that yield a phase diagram with the biologically im- 
portant fluid phase and, additionally, various gel phases. 
Analytical studies have suggested that the phase behav- 
ior of lipid bilayers is dominated by packing effects due 
to the excluded volume of the hydrophobic tails4^ In 
our model, we can draw on the vast knowledge of liquid- 
state theory to control the degree of packing effects and 
local structure of the fluid in order to tailor the weighted 
densities such that the fluid exhibits pronounced packing 
effects. 

The dimensionless, microscopic densities, /5 Q (r), of hy- 
drophilic and hydrophobic species are functions of the 
explicit coordinates of the effective interaction centers 

A*(*) = -jf r )W (9) 

i=l 

where t(i) £ {A, B} denotes the species of bead i. 
The prefactor has been chosen such that the molecu- 
lar density does not depend on the number of interac- 
tions centers per molecule, N. In order to regularize the 
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FIG. 1: The two weighting functions from eqs. and (|12[1 
with a = 0.9AL. The inset shows the derivatives of these 
functions, which are proportional to the non-bonded forces 
between two beads. 

(5-function in the excess free-energy functional of non- 
bonded interactions, Eq. (|2J|, we use a weighted-density 
approximation 4 ^— and define coarse-grained densities 

R 3 nN 

p ma {r) = -^£i«m(|ri-r|)(S Qt(i) (10) 

»=i 

by convoluting the microscopic, molecular density, p a (r) 
with weighting functions, w m . We require that the 
weighting functions are differentiable, vanish for r > AL, 
and are normalized, i.e., J d 3 rw m (|r|) = 1. Liquid-state 
theory for simple liquid a 45 ' 46 as well as integral equa- 
tion theor y 49 ' 50 indicate that it is important to use dif- 
ferent weighting functions to represent the harsh, short- 
ranged repulsion in a liquid and the soft, longer-ranged 
attractions. The second-order terms in Eq. ([2]) typically 
correspond to attractive interactions and the third-order 
terms to repulsions. Therefore, we use different weight- 
ing functions, u>2 and W3, for the second and third-order 
contributions. Both weighting functions are plotted in 
Fig. [T] The longer-ranged weighting function, W2 , con- 
sists of a constant part for r < a and a cubic spline for 
a < r < AL with < a < AL, given by 

{(AL -a) 3 , r < a 
2r 3 - 3 (a + AL)r 2 + 6aALr ■ ■ ■ , 
3aAL 2 + AL 3 , r < AL 

(11) 

A = -15/ [27r(2a 6 - 3a 5 AL + 3aAL 5 - 2AL e )] is a nor- 
malization constant. In the following we use a = 0.9 AL. 
It is used for the mainly attractive, second-order terms. 
The weighting function for the repulsive interactions, W3, 
is the standard choice in Dissipative Particle Dynamics 
models^i 

^ = ^W(AL-r) 2 . (12) 
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It only possesses positive Fourier modes. Negative 
Fourier modes of pair-wise, repulsive interactions give 
rise to cluster-crystallization in dense liquids of soft 
particles^r— Our choice of weighting functions avoids 
the formation of cluster-crystals in the range of parame- 
ters investigated in the following. 

Using Eqs. © and (ITU1) . we rewrite the non-bonded 
interactions in the form 



H 
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which takes the form of a weighted-density 
functional J£ ~ 48 > 55 The density-functional form of 
this coarse-grained interaction free energy controls 
local correlations, e.g., packing effects. Their length 
scale is set by the spatial extent of the non-bonded 
interaction, AL. Unlike density- functional theory, 
however, we obtain the properties not by minimizing the 
density functional but we use density-functional-inspired 
interactions in our soft, coarse-grained model whose 
properties are studied by computer simulation. In this 
way, long-range fluctuations, e.g., undulations of the 
bilayer membrane, are accounted for. 

Finally, we note that weighted densities which give 
rise to strong packing effects deteriorate the quality of 
the mean-field approximation and, consequently, the de- 
coupling between the thermodynamic properties (e.g., 
compressibility and coexistence density) and the liquid 
structure breaks down. Therefore, the model parame- 
ters, p C ocx and K./V, are not identical to the density in the 
hydrophobic interior of the bilayer and its inverse com- 
pressibility. Nevertheless, the approximate equations, Q 
and ([5]), are a useful guide for constructing the model. 



B. Simulation technique 

We applied Multibody Dissipative Particle Dynamics 
(MDPD)££r— to integrate the stochastic equations of mo- 
tion. In MDPD the force Fj acting on each bead, i, con- 
sists of three terms, 



nN 



Here — r,; 



(14) 



r j; and F^ry) = F£(r y ) + F£ b (r y ) 



is the pair-wise, conservative force. The contributions 
from the bonded interactions, F^(ry) are obtained by 
taking the derivative of the potential energy in Eq. (fT]) 
with respect to the coordinates of the beads. 



The non-bonded forces, 



F^ h (r ij), stem from the 



density-dependent Hamiltonian (fT3"|). We rewrite 
Eq. (I13p in a computationally convenient form using the 
expressions for the microscopic and weighted densities. 
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Taking the negative derivative of H n b with respect to , 
we obtain 
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(16) 
(17) 



Ut(i)t(j)a 



Thus, the total non-bonded force is decomposed into a 
sum of pair- wise forces, Fp = ^ . F^. 

The dissipative force, F D (ry , Vy ) , and the random 
force, F R (rij), are used to obtain a canonical ensemble, 
in which the temperature is constant. They have the 



same cutoff, AL, i.e. they vanish for r,j = |ry| > AL. 
For Tij < AL they are given by the DPD form: 



59.60 



F^(r 



F«(r, 



-7W D 0y)(Vy 



(18) 
(19) 



The friction constant, 7, is related to the noise coeffi- 
cient, £, by the fluctuation dissipation theorem, £ 2 = 
2fce7 1 7. 9ij is a stochastic variable with mean, (Oij) = 0, 
and covariance, (9ij(t)9 k i(t')) = {8ij5 k i + SuSj k )6(t - 
t'). The random numbers are drawn from a uniform 
distribution,— and the standard weighting functions for 

DPD 62 



[u,*(r)] 2 =w» 



1-r/AL, r < AL 
0, r > AL 



(20) 



are employed. In the following we use 1 R eo = 3.5 AL. 

Several different thermodynamic ensembles have been 
used in the course of our study. Some simulations have 
been performed in the canonical ensemble (NVT) using 
the standard velocity- Verlet integration scheme with a 
time step of At = 0.005 r.~ Most of the simulations 
have employed an ensemble where the area of the lipid bi- 
layer fluctuated, such that the lateral pressure vanished, 
i.e. Pt = 0. The height of the simulation box L x in 
the direction normal to the bilayer was kept at a fixed 
value. We refer to this thermodynamic ensemble as the 
"NPfT" ensemble. Details of the symplectic integration 
algorithm for this extended ensemble^ are given in Ap- 
pendix [XJ P t — and P s» imply that the bilayer is 
in a state of vanishing mechanical tension, E. The simu- 
lations have been performed by a parallel DPD program 
employing the force-decomposition algorithm devised by 
Plimpton^ 



III. SELF-ASSEMBLY AND BILAYER 
PROPERTIES 

In this section we demonstrate that the lipid molecules 
self- assemble into various morphologies, and we compile 
several static and dynamic properties of the soft, coarse- 
grained model. 
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FIG. 2: Self-assembly of the system p C oox = 17, kN = 
100, xN = 30, N A = 12 hydrophobic (blue), and N B = 4 
hydrophilic beads (yellow). From the initial configuration (a) 
broad wormlike micelles form (b). They coalesce forming a 
bilayer with several pores (c). These pores close slowly and a 
continuous bilayer is formed (d). It happened frequently that 
up to 5 % of the lipids stayed at first in the gas phase, and 
formed after some time one micelle, (a) t = 0, (b) 40 r, (c) 
200 r, (d) 800 r. The red scale bar denotes 2 R co . Created 
with VMD.i££ 



A. Self-assembly 

To study the self-assembly as a function of the molec- 
ular stiffness, we have used three different sets of k s and 
k b (cf. Eq. (fig])); (i) k s = 3.673 and k b = for flexible 
lipids without any bond-angle potential, (ii) k s = 19.0 
and k b = 5 representing lipids with a moderate stiff- 
ness, and (iii) k s = 29.4 and kb — 10 parameterizing 
lipids with a high stiffness. We have explored various 
values of the coarse-grained parameters, 20 < \N < 100, 
50 < kN < 500, and 15 < p C oox < 40, as well as different 
lengths of the hydrophobic tails Na and the hydrophilic 
heads N B with N A + N B = N = 16. We have per- 
formed all simulations in the ATjT-ensemble with £ = 
and have used the same initial configuration comprised 
of n = 1600 lipids in a box of lengths L x = 50 AL, 
and L y (t = 0) = L z (t = 0) = 30 AL. The lipids were 
randomly distributed over the lower half of the box, i.e. 
< x < L x /2 to avoid the formation of multiple bilayers. 

Depending on the parameter set, the lipids self- 
assemble within At < 500 r to one of the following mor- 
phologies: (i) spherical micelles, (ii) cylindrical micelles, 
(iii) wormlike micelles, (iv) bilayers, or (v) bilayers with 
hydrophilic inclusions. A typical pathway of the self- 
assembly of a bilayer is shown in Fig. [5] The results ob- 
tained with other parameter sets are compiled in Tab. HI 



Bilayers form for Na > 11, and inverted structures, i.e., 
bilayers with hydrophilic inclusions, form for Na > 13. 
For kb = 0, 5, 10 we observe only the fluid phase, a 
fluid and a gel phase, and only a gel phase, respectively. 
For \N > 50 wormlikc or cylindrical micelles predom- 
inantly form, whereas for xN < 20 the incompatibility 
between hydrophilic and hydrophobic beads becomes so 
small, that no clear separation between hydrophilic and 
hydrophobic regions is visible. 

The observed sequence of morphologies is consis- 
tent with the geometrical arguments put forward by 
Israelachvili<££ For N A < 11 the amphiphiles have a con- 
ical shape, so that only micelles occur irrespective of 
the other parameters. Na — 11 and Na = 12 result 
in an almost cylindrical shape of the lipids, so that bilay- 
ers form. When the hydrophilic heads decrease in size, 
Na > 12, inverted morphologies appear. With increas- 
ing /o coex each coarse-grained bead interacts with more 
neighbors, so that the mean-field approximation becomes 
more accurate and fluid-like packing effects weaker. This 
marks the crossover to polymeric membranes, where the 
chain number density is typically higher than in lipid bi- 
layer membranes and only a fluid phase is stabled 

The spatial extension of a lipid molecule is of the or- 
der i? eo , but the fluctuations around this mean value are 
largely influenced by kb- The value, kb — 0, corresponds 
to fully flexible molecules and the shape fluctuations are 
of the same order of magnitude as the lipid's size, i.e., 
the conformations resemble a self-avoiding random walk. 
For kb — 10 the lipids are strongly elongated and they 
behave like rods. This gives rise to nematic, liquid crys- 
talline structure of the self-assembled bilayers. 

B. Bilayer properties 

In the following, we focus on the parameter set, Na — 
12, Nb — 4, k s — 19.0, and kb = 5, which gives rise to the 
spontaneous formation of bilayer membranes. The non- 
bonded interactions are set to kJV = 100, %A = 30 and 
Pcoox = 17 or 40. The longest simulation runs lasted 
At w 10 4 t. We used pre-assembled bilayers as ini- 
tial configurations with n = 4680 lipids. In the case of 
Pcoox = 17 two different initial configurations have been 
used - one in the liquid phase, L ai and one in the gel 
phase, Lp. The former one has also been employed as 
initial configuration for p cocx = 40. We compiled the 
obtained properties in Tab. [TTJ 

1. Density profiles 

Stable fluid membranes in a solvent form a bilayer 
structure. The hydrophilic head groups on the outside fa- 
vor contact with the solvent, and the tails constitute the 
bilayer 's hydrophobic interior, which is shielded from the 
solvent. This lamellar structure becomes visible in the 
molecular density profile, which has been recorded sepa- 
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TABLE I: Morphologies observed during self-assembly from a disordered starting configuration 
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^this notation means Na 


= 10, N B 


= 6 










^s: spherical micelles, c: 


cylindrical micelles, w: wormlike micelles, b: bilayer, i: 


bilayer with hydrophilic inclusions 






TABLE II: Static and dynamic properties for kN 


= 100, X N = 


30, n = 4680 












Pcoex = 40 


pcoex — 17 


Pcoex — 17 










Fluid (L a ) 


Fluid (L a ) 


Gel (Lp) 


Area 




(A) [RL] 


63.7(1) 


109.6(1) 


99.7(1) 


Area Compressibility 




k A [lQ- A Rljk B T] 


4.17(1) 


3.53(1) 


0.14(1) 


Area per Lipid 




(« 


) [W~ 2 R 2 co ] 


2.722 


4.684 


4.262 


Bulk Density 




PA [l/Rl] 


42.9(1) 


22.4(1) 


23.6(1) 


Width of Hydrophobic Layer 


W 


[Rco] 


1.21(1) 


1.42(1) 


1.53(1) 


Total Thickness 




t 


[Reo] 


1.63(1) 


1.87(1) 


2.05(1) 


Aspect Ratio 




Vj/JQj [1] 


7.3 


6.5 


7.4 


Bending Rigidity (spect 


rum) 


K 


[k B T] 


19(1) 


15(1) 




Bending Rigidity (from 


k A ) 


K 


[k B T\ 


13 


21 




Molecular Diffusion Constant 


D 


[RlJr] 


1.5(1) • 10~ 3 


5.0(1) ■ 10" 4 


2.6(1) ■ 10" 7 



rately for the two leaflets. To avoid a broadening of these 
profiles by thermal undulations, the bilayer has been sub- 
divided laterally into small cells of size, AL x AL. In 
each cell, the local bilayer position has been determined 
and the profiles have been averaged with respect to this 
local bilayer position over all cells and along the trajec- 
tory. As noted above, packing effects result in a deviation 
of the density of ^4-beads in the hydrophobic core from 
the estimate of the coexistence density, /3 C oex, provided 
by mean- field theory. The width, w, of the hydropho- 
bic core is estimated in both phases by measuring the 
full width at half maximum (FWHM) . The total bilayer 
thickness, t, is the distance from the center of mass of 
the hydrophilic head groups on the one side to that of 
the apposing side. 

Fig. [3] depicts two density profiles, p a (x), across the bi- 
layer. They correspond to the liquid and gel phases, L a 
and Lp, which have been observed at /9 C oox = 17. Both 
show separated peaks for the hydrophilic heads and the 
hydrophobic tails, which implies that the coarse-grained 
lipids indeed form bilayer membranes. A closer inspec- 
tion of the densities of the hydrophobic interior shows, 
that the two leaflets are clearly distinguishable, but that 



there is no dip in the center of the density profile, like it 
is known from atomistic or systematically coarse-grained 
models including solvent. In fact, we find a flat profile in 
the fluid phase and a hump in the gel phase, the latter 
being caused by an overlap of the last bead of the lipids 
from each side (not shown). 

The remaining overlap between the apposing leaflets in 
the liquid phase and the hump in the density profile in the 
gel state might arise from three different reasons: (i) The 
molecular shape in our model is rather finely discretized 
and the lipid tails are rather flexible. If the molecu- 
lar shape becomes more rod-like, the density profile at 
the center is expected to develop a dip due to molec- 
ular packing. This could be achieved by a decrease of 
the number of beads per lipid or an increase of the bond 
stiffness, (ii) If the incompatibility between hydrophobic 
and hydrophilic segments increases, the bilayer thickness 
will increase and the intcrdigitation between the appos- 
ing leaflets will decrease, (iii) The flat density profile 
could also arise from the lack of solvent molecules. Since 
there is no solvent exerting pressure on the membrane, 
the lipids might have to interdigitate slightly, so that the 
whole bilayer remains stable. This would indicate a gen- 
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FIG. 3: Molecular density profiles p a (x) across the bilayer 
with respect to the local midplane at p COC x = 17. The upper 
graph shows the density profile of the Lp phase, while the 
lower one presents the L a phase. The dotted lines denote the 
densities of the hydrophilic beads, the dashed lines indicate 
the densities of the hydrophobic beads for each leaflet, and 
the solid line marks the sum of the two hydrophobic densities. 
Finally, the dash-dotted line indicates p COC x, which has been 
used to derive the interactions of the hydrophobic beads. 



eral problem in solvent-free models. We are not aware, 
however, of density profiles in the gel phase for another 
solvent-free model. 

Although the width, w, of the hydrophobic core and 
the area per lipid, (a) = (A)/(n/2), depend on the de- 
tails of the soft, coarse-grained model, the dimension- 
less aspect ratio, w/ •*/ (a), can straightforwardly be com- 
pared to experiments. The most common two-tailed 
lipids have aspect ratios in the range of w/ \J (a) ~ 3 — 5, 
whereas our simulations for single-tailed lipids yields 
6 < w/y/(a) < 7. If one assumes that two single-tailed 
lipids of our coarse-grained model glued together result 
in one two-tailed lipid, one will double the mean area per 
lipid, (a), and obtain an additional factor of \[2 in the 
denominator of the aspect ratio, thereby, obtaining as- 
pect ratios that are in good agreement with experimental 
values. Noteworthy, it is impossible to obtain the right 
aspect ratio simply by selecting different interaction co- 
efficients at fixed discretization, N, and fixed molecular 
architecture. At this level of coarse-graining the latter is 
important and must be taken into account, if one tries to 
map a specific kind of lipid. In this study, however, we are 
content with the simplest molecular architecture, i.e. lin- 
ear molecules. We establish a conversion factor between 
the unit of length in the simulation and in experiments 
using the mean area per molecule. For synthetic lipids 
like DPPC, DPPE, or DLPC (a) « 60 A 2 £L™. so that 
we find in the fluid phase at p cocx = 17 an equivalence of 
1 R eo = 3.8 nm. 



2. Elastic properties 

In the TVPtT-ensemble the projected area of the bi- 
layer, A, is fluctuating. These fluctuations are related to 
the area compressibility, ]$a , vi a 15 ' 69 ! 71 



k A =p 



(a 2 ) - (Ay 

(A) 



(21) 



k- 1 

A ' 



This formula neglects undulations of the bilayer, which 
result in a difference between the projected and the true 
surface area of the bilayer^ 9 - However, for the small 
patches of a bilayer used in this study, this difference 
is negligible. 

We have measured Ua according to Eq. f[2T|) for the 
fluid and the gel phase at p CO cx = 17 as well as the fluid 
phase at p cocx = 40. Using the conversion factor for 
the unit of length from above and converting the area 
compressibilities to area compression moduli, Ka = 
we find Ka = 79.4 mN/m for the fluid phase at p c 
17, 67.2 mN/m at p cocx = 40, and 2045 mN/m in the 
gel phase at p cocx — 17. Typical experimental results for 
two-tailed lipids in the fluid phase yield values of Ka ~ 
240 mN/m£ We attribute the lar ger area fluctuations 
of our model bilayer to the softer interactions and the 
reduced number of degrees of freedom. 

Another important quantity is the bending rigid- 
ity, k, which measures the cost of undulations of the 
bilayerj 7 ^— It is frequently calculated in particle-based 
simulations. 14 ' 15 ' 35 ' 71 i 76 ~ — If the fluctuations are small, 
the free energy of a curved membrane is given by the 
Helfrich-Hamiltonian 



H 



d z r 



n{V 2 hf +E(V/i) 5 



(22) 



where h(r) in the Monge gauge is the bilayer's height 
above a reference plane. By inserting the Fourier expan- 
sion, h(r) — /i q exp (iqr) with q = 2ir(n x , n y )/L, in 
Eq. (f2"2")) . we observe that the different modes, /i q , decou- 
ple, and, using the equipartition theorem, we obtain for 
the power spectrum 7 ^ 



A( K q 4 + Hq 2 )' 



(23) 



Eq. ([2"3")) has been derived for the canonical ensemble, 
i.e., a constant box size. Here, we also used it to fit our 
data in the A^P t T-ensemble. Since the lateral box lengths 



fluctuate in this ensemble, (\h c 



denotes the mean in- 



tensity of each mode at the time-averaged wave vector 
q = n • (2ir/L{t)). Additionally, we included the spectral 
damping factor in the calculation of (|/i q | 2 ) that arises 
from the interpolation of the continuous bilayer position 
onto a grid; 14 From a computational point of view, we 
calculate the bilayer position, h(r), on a quadratic 16 x 16 
grid by averaging over the perpendicular distances of all 
hydrophobic beads from the reference plane. Then we 
perform an FFT of h(r) yielding the amplitudes /i q . 
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perimental units. Using 1 R co = 3.8 nm and a typ- 
ical lipid diffusion coefficient at room temperature of 
D = 5 /j,m 2 /sf£ we obtain 1 r = 1.5 ns. It is interesting 
to relate this identification of time scale in our coarse- 
grained model to the occurrence of flipflop events. Unfor- 
tunately, the flipflop rate hardly deviated from zero; we 
have observed only a very small number of events even 
in the longest simulation runs. Therefore only a lower 
bound for the mean time, (t), between two flipflop events 
is presented here. We find (t) > 10 8 r « 0.15 s. This 
is reasonable, since a passive flipflop is a thermally ac- 
tivated process, which happens on an experimental time 
scale of 1 event per molecule per day 84 i 85 

IV. MAIN PHASE TRANSITION 



FIG. 4: Power spectrum of the height fluctuations (\hq\ ) of 

the fluid phase at p COGX 17. Inset: k B T / A(\hq\ 2 )q 2 plotted 

as a function of q 2 . For small q 2 the data points are fitted to 
a straight line through the origin with slope k. 



Fig. [4] shows the power spectrum of bilayer fluctua- 
tions in the fluid phase at p cocx = 17 in a tensionless 
state. By fitting kBT/A(\hq\ 2 )q 2 as a function of q 2 to 
a straight line through the origin, we extract k from the 
slope (cf. Eq. ([23]) ). Experimental values of the bend- 
ing rigidity k lie for most biological membranes within a 
range of 5 — 60 ksTj^ and our results match this order of 
magnitude for both fluid systems under study. However, 
no bending rigidity could be obtained by this method in 
the gel phase. 

Another independent, but rather crude estimate of k 
is provided by the area compressibility. It has been sug- 
gested thafci&^S.Zi^ 



b-k A ' 



(24) 



where t is the thickness of the bilayer. There has been 
some debate about the value of the geometric factor b. 
Here we use b = 48. The obtained values match the order 
of magnitude of the values extracted from the undulation 
spectra. 



3. Diffusion 

Finally, we have measured the lateral mean-square dis- 
placements of the lipids' center of mass, and obtained the 
two-dimensional self-diffusion coefficient, D, from 



D = lim 

t^oo it 



-/(r CI 

It V 1 



it) 



(25) 



This equation neglects all undulations but it is a good 
approximation for the systems under study 81 ' 82 

We used D (see Tab. |H|) to establish a mapping be- 
tween the unit of time, r, in the simulation and ex- 



Depending on the control parameters, nN and p COC x, 
the lipids self-assembled into bilayers of different ther- 
modynamic phases. We observe the fluid phase, L a , the 
non-interdigitated gel phase, Lp, the fully-interdigitated 
gel phase, Lpi, and a tilted gel phase, Lp/M Among 
the different phase transitions, the main phase transi- 
tion, Lp f-> L a , is definitely the most-important one. It 
has many of the characteristics well known from first- 
order transitions, like pronounced hysteresis effects, the 
occurrence of metastable states, and sharp peaks in the 
response functions. 

We used three different, but not independent methods 
to locate phase coexistence. First, we have applied a com- 
bination of Umbrella Sampling (US) and the Weighted 
Histogram Analysis Method (WHAM) to compute the 
free energy ) 23 i 24 ' 27 ~ — F(p coex ), in the vicinity of the main 
phase transition. Second, we have utilized Free Energy 
Perturbation theory (FEP) to extrapolate the free en- 
ergy branches of each phased Finally, we have used a 
histogram reweighting scheme to calculate the specific 
heat, C(p coex )£iiZ. 



A. Order parameters 

Several order parameters characterize the main phase 
transitions ; 14 ' 16 ' 86 We chiefly employ the orientational or- 
der parameter 




(26) 



where cos fljj+i = n- (r J+ i — r 3 -)/ |r J+ i — r^] denotes the 
angle between the local, normal vector, n, to the bilayer 
and the bond vector between two succeeding hydrophobic 
beads j and j + i sums over all molecules in the bilayer 
and the average is taken over an ensemble of bilayers. 
5=1 means that all lipids are perfectly aligned parallel 
to n, S — indicates isotropically distributed directions, 
and in the case S = — 1 the lipids are perfectly aligned 
in the plane of the bilayer. Since n(r,) is a function of 
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the coordinates of many lipids, its calculation involved a 
triangulation procedure^ where we described the bilayer 
midplane by a set of small triangles with a unique normal, 
n(rj), in each triangle. 

In the gel phase the lipids form a two-dimensional 
structure with 6-fold symmetry. We probe this inter- 
molecular packing by the order parameter 




Here, rii denotes the number of lipids adjacent to lipid i 
(as determined by a Voronoi tesselation) , and 4nj the an- 
gle between the vector from the center of mass of lipid i to 
that of lipid j, and some arbitrary but fixed direction in 
the plane of the bilayer. ipQ = 1 indicates perfect hexag- 
onal symmetry over the entire bilayer, whereas ^6 = 
signals the absence of bond-orientational order. 

Both order parameters, S and tp6> clearly distinguish 
between the liquid and the gel phase, but they differ 
in one crucial point: S is composed of additive contri- 
butions, which only stem from conformational, single- 
molecule properties and therefore its change in response 
to moving a segment can be easily computed. The oppo- 
site is true for ipg, which only reflects bond-orientational 
order caused by intermolecular packing and requires 
the computationally intense Voronoi tessellation in or- 
der to identify the neighbors of a lipid. The main 
phase transition simultaneously involves both, a change 
in the in-plane degrees of freedom that dictate the bond- 
orientational order and a change in the conformational 
degrees of freedom. 44 

We choose S as the single reaction coordinate (or- 
der parameter) for the liquid-gel transition, because the 
conformational and the bond-orientational transition are 
coupled, and S is considerably easier to compute than 

-06- 



B. Determination of coexistence point 

Here, FEP has been used to calculate the free-energy 
difference between two systems that only differ in their 
non-bonded interactions, or more precisely, that only dif- 
fered in the parameter, p cocx . To this end, we sample 
configurations at a reference density, po, and calculate 
the free energy difference, F(p) — F(po), with respect to 
a system with a different density, p, by 2 ^ 

F(p) - F{ Po ) = -k B T ln(e X p(- /3 A-Hnh{p))) o - (28) 

Here (• • • ) stands for an ensemble average of the refer- 
ence system and AH n b(p) is the difference of the non- 



bonded energies between these two systems, 

AH n b(p) = H n b{p) ~ 'Hnb(Po) 
= [Vafl(p) ~ V a p(po)} 

+ [w a/3l (p) - w Q/ 3 7 (po)] (29) 

where the integrated densities, P a p and Q a f3y are defined 
by 

P a /3 = ^2s at ^)p 2 (r l ) (30) 

i 

Q a p-f = 8 a t(i)fap h-y i r i) ■ ( 31 ) 

i 

Since FEP samples only the phase space of one thermo- 
dynamic phase, we cannot locate the coexistence of two 
phases. However, it is well suited to explore the free en- 
ergy branch F{p) of a single phase. 

In contrast to FEP, the combination of US and WHAM 
allows a direct location of the phase coexistence. At first 
the free energy profile, F(S), as a function of the order 
parameter, S, is calculated for a specific set of expansion 
coefficients. Let Sn and S se \ denote the order parameter 
in the fluid and in the gel phase, respectively. To obtain 
F(S), bilayer configurations have to be uniformly sam- 
pled for all values of S in the interval Sg < S < S ge \. 
However, the unfavorable configurations in the misci- 
bility gap are unreachable by conventional Boltzmann 
sampling because their statistical weight is exponentially 
small. By including an additional US potential, Wi, we 
force the system to sample also these unfavorable config- 
urations. Specifically, we add the harmonic potential 

W i = ^n(N A -l)[S-S i } 2 . (32) 

that biases the simulation to keep S in the vicinity of 
St. Here k — 1 k B T is a spring constant that measures 
how strong deviations from Si are penalized. We have 
used an equidistant spacing of the Si with ASi = 0.01 in 
the range 0.15 < Si < 0.8 to sample the whole interval 
uniformly. 

A simulation has been performed for each value of Si, 
in which we have recorded a trajectory of the order pa- 
rameter, S, the total energy, U, and the integrated densi- 
ties, P a /3 and Q a p-y These quantities are used to reweight 
the trajectories to different values of p CO cx (cf. Eq. [29]) . 
Each of them is binned into a normalized histogram that 
measures the biased probability density of visiting S in 
a run with potential, Wj. In the subsequent weighted 
histogram calculation, this bias is removed from the his- 
tograms and all individual histograms are combined into 
one unbiased histogram, po(S), in a way that the statis- 
tical error is minimal. For brevity, we omit the computa- 
tional details and refer to the original work . 27 ' 28 Once the 
Boltzmann probability distribution, Po(S), is available, 
the free energy, F(S)/ k B T, is computed as the negative 
logarithm. 
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FIG. 5: Hysteresis loops of the order parameters S and ^6 
for kN = 50 . . . 125, \N = 30. At each step, the bilayer was 
simulated for At = 100 r. The step size between two simula- 
tions was Ap coex = 1.0, and the arrows mark the direction, in 
which pcoex was proceeded. The dashed line at p COC x = 17.27 
indicates the transition point for kN = 100. 



H 
M 



O.i 



0.4 



0.0 



-0.4 



0.2 



0.1 



0.0 



1 1 1 1 1 1 

Fluid Phase L 

_ a 






A 


- 




v AF=648 k B T - 


- 


p = 17 


N. Gel Phase 

v \_L—^ 




I , I 


i 




0.2 0.4 

: S 







— Pure Fluid 
Pure Gel 

— Reweighted 



16.8 



17.0 



17.2 



17.4 



FIG. 6: Inset: Free energy, F(S), at p CO ex = 17. There was 
a difference of AF = 648 fcsT between the minima of the L a 
phase and the Lp phase. Main Panel: F(p C oex) obtained from 
histogram reweighting in comparison to two independent FEP 
calculations of a pure L a and a pure Lp phase, whose offset 
AF at pcocx = 17 was known from the inset. Both curves 
intersected at pfep = 17.26 (dashed gray line). 



F(p, S) has been computed by histogram reweighting 
similar to Eq. (|28p . In contrast to the FEP calculations of 
the pure phases, F(p, S) includes contributions from both 
phases. By taking the integral of F(p, S) over all S we 
obtained F(p); a phase transition in this quantity is vis- 
ible as a point with a rapidly varying derivative because 
finite size effects lead to a rounding of the transition. 

We have also applied the reweighting procedure to U, 
so that the probability distribution p(p, U) becomes avail- 
able. From this quantity we calculate the mean total 
energy, (U)(p), and the specific heat 

c ((*) - W) 

hs ~ (k B TT • [ ' 

which serve to locate a first-order transition. 

A first estimate of the position of the main phase tran- 
sition is obtained from the center and the width of a hys- 
teresis loop. Therefore we have simulated pre-assembled 
bilayers with 1600 lipids that were initially in the fluid 
phase at p C ocx = 40. We have performed several succeed- 
ing cycles with p coex running from 1 1 to 40 and vice versa 
in steps of Ap cocx = 0.5 or 1.0 for kN = 50, 75, 100, 125. 

Near the main phase transition, large hysteresis effects 
occur in S and ipQ as shown in Fig. [5] The loops for 
both order parameters differ only quantitatively. They 
are weakly shifted, and intramolecular order persists up 
to slightly higher p coex than the intermolecular order. 
The widths of the loops grows with increasing kN , i.e., 
metastable domains persist up to higher p CO cx. Addi- 
tionally, the amplitudes of the order parameters increase 
indicating different thermodynamic phases. For instance, 
in the case kN = 125, two distinct gel phases {La and 
Lgi) occur. Their transition is visible as a dip in both 
order parameters near p coex ~ 20. 



We focus on the system, kN — 100 and \N = 30, 
where the results in Fig. [5] have indicated that the main 
phase transition is located in the interval 15 < p CO cx < 
18. We calculate F(S) by means of US/WHAM from 
simulations at different coexistence densities, p CO cx = 
16.28, 17.00, 17.29, with 2-5 different initial configu- 
rations in the 7VP t T-ensemble for At = 3000 t. It is 
advantageous to start the simulation from an initial con- 
figuration where both phases are already present^ 

The inset of Fig.[S]shows F(S) at p C oox = 17. The two 
visible minima correspond to the metastable L a and the 
stable Lp phase. However, the offset AF = 648 /c.bT be- 
tween these minima indicates that the gel phase is ther- 
modynamically stable. To locate the phase transition, 
we reweight p C ocx searching for a rapid variation (i.e., 
rounded discontinuity) of the slope that signals the phase 
transition. Such a kink occurs at pus — 17.29 (cf. main 
panel of Fig. [5]) indicating the crossing of the free energy 
branches of the different phases. 

It is convenient to introduce a normalized order pa- 
rameter, AS = (S — Sg)/(S se \ — 5fl), so that the minima 
of the free energy in the fluid phase at Sn and in the 
gel phase at 5 ge i correspond to AS — and AS = 1, 
respectively. Fig. [7] depicts F(AS) at p COC x = 17.29. At 
this point both phases have equal statistical weight, and 
they are separated by a free energy barrier with a plateau 
value of F s iab = 151.0(5) fcgT. We calculate AS from the 
abscissae of the minima, Sr — 0.212 and S gc \ = 0.625. 

To confirm the transition point, we employ Eq. (|33[) 
to calculate C(p coex ) and (Z7)(p coex ) by reweighting (see 
Fig. [8]) . The main phase transition is visible as a sharp 
peak in C(p coeK ) at psh = 17.27, as well as a steep rise 
in (C/)(p COC x) at the same position. The slow rise of (U) 
in the interval 16.8 < p CO cx < 17.27 can be attributed to 
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FIG. 7: F(AS) at phase coexistence compared to the phe- 
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FIG. 8: Specific heat C (above) and total energy (U) (below) 
obtained from histogram reweighting as a function of /9 C ocx- 
The main phase transition is visible as a sharp rise in ([/), 
but also as a peak in C at psH = 17.27 (dashed red line). The 
dotted gray lines mark the points where simulations using US 
have been performed. 



the gradual melting of the hydrophobic tails. 

Finally we have conducted two additional, indepen- 
dent simulations at p CO cx = 17 in the iVPfT-ensemble 
without an US potential. One initial configuration was 
prepared in a pure fluid phase (L ai metastable) and 
the other was prepared in a pure gel phase {Lp). The 
free energy branches of each phase are extrapolated with 
FEP (cf. Eq. (|28]1). In this method the relative free 
energy difference between both branches remains unde- 
termined. However, it has already been computed by 
the offset between the branches at p C oox = 17 yielding 
AF = 648 k B T. The main panel of Fig. [5] depicts the 
two, correspondingly shifted branches of F(p cocx ), which 
intersected at pfep = 17.26. 



Gratifyingly the US /WHAM results for the free energy 
F(AS) and F(p cocx ) are consistent, indicating the high 
statistical accuracy of our data. In Fig. [Jo] we present 
the two branches of the free energy, F(p cocx ), of each 
phase obtained from FEP in comparison to the result 
from US/ WHAM. In the fluid phase both methods com- 
pletely agree, however, in the gel phase there is a small 
difference discernable. Between p coe x = 17.15 and 17.29 
the FEP calculation slightly overestimates the free en- 
ergy by AF = 0.02 k B T per lipid, which arise from a 
gradual loss of bond-orientational order as the transition 
is approached from the gel phase. Therefore the result, 
Pfep, is less accurate than the other estimates. 

A similar way of determining the phase coexistence 
point has been applied earlier J^- In that study the rela- 
tionship between the branches is fixed by the knowledge 
of the two bulk free energies, which were extracted from 
mean-field theory. A similar calculation is possible in the 
fluid phase of our coarse-grained model, but it is not ac- 
curate in the gel phase, where correlations between the 
lipids are essential. These correlations, which are cap- 
tured by in our simulations, are clearly visible, e.g., in 
the order-parameter tpQ. 

The three estimates of the location of phase coexis- 
tence, pus, pfep, and psh, nicely agree with each other. 
The main error source of our estimate of the phase co- 
existence, however, stems from possible sampling error 
along the US path that reversibly connects the liquid and 
the gel phase, which are difficult to estimate. The consis- 
tency of the results suggests that the liquid-gel transition 
for the parameters, kN = 100 and \N = 30 , occurs at 

Pc * ocx = 17.3(1). (34) 

Thus, the uncertainty in the location of the transition 
point is reduced by a factor of 50 compared to the un- 
certainty in the hysteresis loops. 

C. Bilayer configurations 

Besides the point of the phase coexistence, F(AS) 
also offers an insight how typical configurations of the 
finite bilayer inside the miscibility gap look like^Zr— Let 
us consider a lipid bilayer in the fluid phase. If A5* is 
slightly increased from the value AS* = this small in- 
crease will be distributed homogeneously throughout the 
bilayer. The excess free energy of this undercooled fluid 
bilayer up to second order in AS is given by a Taylor 
expansion around the minimum 

F uf = ^ (AS) 2 , (35) 

where fc u f is a constant measuring the response of the 
system to changes in AS". 

In a macroscopic system an undercooled bilayer 
is metastable and the lipids will condense into two- 
dimensional droplets of radius, R, that consist of the 
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FIG. 9: Typical bilayer configurations inside the miscibility gap (duplicated across each periodic boundary), (a) shows gel 
droplets in a fluid phase (AS = 0.26), (b) shows the slab geometry (AS = 0.60), and (c) shows fluid droplets in a gel phase 
(AS — 0.84). Each lipid is colored by its local order parameter: the Lp phase is white, the L a phase is dark, and intermediate 
values are gray. Created with VMD.— (d-f): Voronoi tesselation^ of the same configurations as before, showing the area 
occupied by each lipid. The lipids in the Lp phase ordered in a hexagonal structure, while no such structure is present in the 
L a phase. 



thermodynamically stable, gel phase (cf. Fig. c). In 
the framework of classical nucleation theory, the excess 
free energy of such a droplet is given by the droplet 's 
perimeter and the thermodynamic line tension, a, i.e. 

F drop = 2tt<jR. (36) 

Note, however, that the thermodynamic line tension de- 
pends on the length scale, i.e., the perimeter of the drop. 
Since the fluid and the gel phase have both the same 
free energy at coexistence, there is no bulk contribution 
to Eq. (|36| from the interior of the droplet. Since the 
lipids occupy in both phases roughly the same area (cf. 
Tab, mi), the area of the droplet is in good approximation 
proportional to AS, i.e. 

irR 2 ~ L 2 AS. (37) 

where we have used that the normalized order-parameter, 
AS, quantifies the fractional area of the gel phase. Com- 
bining Eqs. fl35) and (|37]). one obtains R ~ AS 1/2 and 

i^drop = 2<7LVttAS. (38) 



If the two-dimensional droplet grows larger, its size will 
become comparable to the linear dimension, L, of the 
simulation box. Then it is more favorable to form a gel 
phase slab that is separated from the fluid phase by two 
plane interfaces of length L (cf. Fig. [5p). In this case, 
the excess free energy is independent of AS, i.e. 

F Blab - 2aL (39) 

Increasing AS even further, one observes the reverse set 
of configurations. The slab of the fluid phase grows thin- 
ner and thinner, and at some point it becomes favorable 
to form a fluid droplet surrounded by the gel phase. The 
radius of the fluid droplet decreases while AS increases. 
Finally, the droplet vanishes and the lipids form a homo- 
geneous gel phase. 

A fit of F(AS) to Eq. (|35l) in the vicinity of the minima 
yields the two constants, A; u f = 6750 ksT and k og = 
4290 fcflT, that quantify the response of the bulk phases 
to changes in AS. We indicate the resulting parabolas 
in Fig. [71 In addition, we also plot AF slab from Eq. (j3"5j) . 
as well as Eq. (|38p for the droplet shape on each side of 
the free energy profile. 
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FIG. 10: For AS < 0.5 only a small amount of the lipids 
had a straight conformation, and it was unfavorable to form 
gel domains; no hexagonal symmetry, measured by tpe, was 
visible. For AS > 0.5 enough lipids acquired such a confor- 
mation, so that gel domains formed and hexagonal symmetry 
set in quickly. 



FIG. 11: To calculate the Gibbs dividing surface the bulk val- 
ues of the order parameter Sgei and Sa are computed (dashed 
lines). After that, the integral criterion is applied in a narrow 
interval h — Sy ■ ■ ■ h + Sy surrounding each edge of the slab. 
The position of the interface h is chosen, such that Af — A g . 



V. LINE TENSION 



It is interesting to note that the left half of Fig. [7] 
with AS < 0.5 fits the phenomenological expressions, 
Eqs. (J2S1) an d well, while the data for higher val- 

ues of the order parameter exhibit larger deviations. 
In addition, fc u f > fc og , i.e., the gel phase, Lp, has a 
smaller response with respect to changes in AS" than the 
fluid phase, L a . This discrepancy stems from the onset 
of bond-orientational order. Fig. [10] shows that bond- 
orientational order is weak for AS < 0.4. For larger val- 
ues of the order parameter, however, the lipids acquired 
also bond-orientational order, which results in a decrease 
of F(AS). The minimal free energy is finally reached 
in a state with hexagonal symmetry of the lipids. Only 
the trailing end beads of the lipids interdigitate with the 
ones from the apposing leaflet and constitute a thin, dis- 
ordered layer at the center of the bilayer. 

Additional deviations from the simple phenomenologi- 
cal estimates arise from the interaction between the lines 
that separate the liquid and the gel domains. For in- 
stance, the thermal fluctuations of the two lines in the 
slab geometry induce attractive Casimir forces^ While 
the slab is growing thinner and thinner at the crossover to 
the droplet geometry, these Casimir forces become more 
pronounced and reduce the free energy. However, there 
is an additional reduction of the free energy coming from 
transversal line fluctuations^ These become pronounced 
for small widths of the slab and finally lead to its de- 
struction. While it is in principle possible to study these 
fluctuation mediated interactions by careful inspection of 
F(AS) at the edges of the plateau, we did not investigate 
these effects in further detail. 



In this section we present two different methods for 
obtaining the line tension from the slab-configurations in 
the middle of the miscibility gap, extracting the bare line 
tension, A, and the thermodynamic line tension, a. 

A. Bare line tension, A 

If the boundary line separating the two domains in 
the slab geometry is smooth and free of overhangs, one 
can describe its position by a function, h(z). The sta- 
tistical properties of h(z) follow from the capillary wave 
Hamiltonia n 94 ' 95 

l z 

H cap = AL 2 + ^ Jdz [h'(z)} 2 , (40) 
o 

where L z is the projected length of the line. 

Routinely, A is determined from h(z). To this end, 
one expands h(z) in a Fourier series and investigates the 
power spectrum of fluctuations, (\h q \ 2 ). A fit to the ex- 
pression 



which one derives from Eq. (|40p using the equipartition 
theorem (cf. Appendix IB"]) . then yields A. 

Several different schemes to locate the position of 
such an interface are knownj2& We use an integral 
criterion ) 97 ' 98 in which we subdivide the profile of the 
local order parameter, S(y,z), into N z = 16 horizontal 
stripes with a width of Az = L z /N z . Each stripe is 
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binned into two histograms, one for each leaflet, with a 
bin width of Ay = 2 AL (where AL denotes the range 
of the non-bonded interactions), yielding in total 2N Z 
histograms per snapshot. The bulk values of the order 
parameter in the gel phase, S se \, and in the fluid phase, 
<Sfl, have been extracted once for each snapshot. As il- 
lustrated in Fig.[TTJ h(zi) in stripe i has been calculated 
for each side of the slab separately as the position of the 
Gibbs dividing surface, such that 

h h+Sy 

J dy (S(y, zi) -S a )= J dy (S ge} - S(y, *)) . (42) 

h—5y h 

Here Sy = 10 AL is used to confine the integral to a 
narrow region surrounding the Gibbs dividing surface, 
so that fluctuations of the bulk influence the position 
of the interface minimally. Once the function, h(z), is 
computed, the fluctuation power spectrum (|/i q | 2 ) is cal- 
culated via FFT and averaged. 

We compared two different ways of calculating this av- 
erage. On the one hand, we argue that the bilayer is 
essentially a two-dimensional object, so that only the 
boundary lines at different sides of the slab, but not 
on different leaflets, fluctuate independently. This leads 
to an average calculated over two independent line con- 
figurations per snapshot. On the other hand, we can 
also locate the boundary in each leaflet independently 
and study their fluctuation spectra. Hence, the aver- 
age included four different lines per snapshot. Finally, 
each average was divided by the spectral damping fac- 
tor [sin(7rn/iV z )/(im/N z )] 2 J^- For small lateral wave vec- 
tors, we expect that the boundaries in the two apposing 
leaflets are coupled and both methods yield the same, 
bare line tension in the limit, q — > 0. To record the 
fluctuation spectrum of the interface line a rectangular 
shape of the simulation box is chosen, with L y > L z . In 
this asymmetric situation, the slab will attain the lowest 
intcrfacial free energy if it aligns parallel to the z-axis. 
So the orientation of the interfaces is dictated by the 
system geometry. Specifically, we have assembled an ini- 
tial configuration with n — 14040 lipids, at p COCK = 17 
having L y — 3L Z w 30 R eo . Note that the length of 
the interface, L z , is the same as in section HVl only the 
other box length, L yi has changed. This configuration 
is simulated in the ./VPtT-ensemble with two indepen- 
dent degrees of freedom. L y and L z . allowing anisotropic 
area fluctuations. Additionally, we employ a US poten- 
tial with So = 0.34, k — 5.0 fc^T, driving the system to 
the desired slab geometry. After an equilibration time of 
At = 10 4 r the box lengths fluctuate around mean values 
of (L y ) — 31.78 R eo and (L z ) = 10.0 R eo , and both inter- 
faces on both leaflets of the bilayer are flat. Unlike the 
experimental situation, where only domains of spherical 
shape are observed, there is no Laplace pressure because 
the interfaces are not curved and thus the two phases 
coexist at the same vanishing lateral pressure. 

Fig. [T3] shows the line fluctuation spectra for the two 
different ways of averaging with two and four indepen- 
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FIG. 12: Power spectrum (|/i<j| 2 ) of the boundary line fluc- 
tuations calculated for 2 (solid) and 4 (dotted) independent 
interface lines with the standard mean error in comparison 
to the value calculated from F(AS) (dash-dotted). The 
dashed grey lines indicate the estimated value of the UV- 
cutoff g m ax = 2n / a. Inset: (\h q \ 2 )q 2 plotted as a function 
of qL z . For small q this function became constant and inter- 
sected the y-axis at ksT/\L z . 

dent lines. In the inset of Fig. [12l (\h q \ 2 )q 2 is plotted as 
function of (qL) 2 . For small q this expression becomes 
linear and intercepts the y-axis at ksT / XL Z . Fits yield 
A 2 = 5.17 kBT/R eo for two lines and A4 = 4.35 kBT/R co 
for four lines, respectively. 

For small wave vectors, q, both graphs in the fluctua- 
tion power spectrum of Fig. [12] asymptotically approach 
the same q~ 2 power-law. This indicated that interface 
fluctuations are correlated on large scales. At larger val- 
ues of q, the lines fluctuate independently and, concomi- 
tantly, the high-g estimate of the line tension is lower 
for the data extracted from four lines than the data for 
two lines. For qL z /2n > 4 the graph shows clear de- 
viations from the simple power-law, which is expected, 
because the description of the line by the capillary wave 
Hamiltonian breaks down on microscopic scales. Thus, 
L z s» 10 R co is too small to observe the expected q~ 2 - 
scaling over an extended g-range. We compute the final 
result by taking the average of both values A2 and A4 
and the deviation as the error bar. This yields the final 
estimate of the bare line tension, A = 4.8(4) /csT/i? eo . 



B. Thermodynamic line tension, a 

A different measure of the free energy cost of a phase 
boundary - the thermodynamic line tension, a - is ex- 
tracted from the free energy profile F(AS) in Fig. [7] The 
excess free energy of the slab configuration is dominated 
by the interfacial free energ y 87 ' 99 , F s i a b = 2aL z , provided 
that the system is large enough, for the two interfaces not 
to interact. _F s i a b is readily obtained from the plateau 
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value of F in the center of the free energy profile. How- 
ever, there is a difference between the bare line tension, 
A, and the thermodynamic line tension, a, which depends 
on the length scale. It is shown in Appendix IBl that the 
two quantities are related via 

a = X + n max In (2nf3XL z ) + 2 In (n max !) , (43) 

which describes the renormalization of the bare line ten- 
sion A by fluctuations. n max denotes the closest integer to 
L z /a and a is a UV-cutoff that characterizes the smallest 
scale, on which the fluctuations of the line are describable 
by a capillary-wave Hamiltonian. a has to be indepen- 
dently determined. 

To calculate the bare line tension, A, from the height 
of the free energy barrier, we estimate the value of the 
UV-cutoff a graphically from the intersect of the q~ 2 
power-law at small q and a constant fluctuation strength 
at higher q. We find g max = 2n/a = 1.69 R^o , i.e., 
a ~ 3.7 i? 00 and thus n max = 3. By taking the 
plateau value F s i a b = 151.0(5) fcgT from Fig.[7]we obtain 
a = 7.55(3) kBT/R eo , and by solving Eq. (|43|) numeri- 
cally for A, we obtain Xf = 5.45(3) fcsT/i? eo , which is 
in good agreement with the estimate obtained from the 
spectrum, A = 4.8(4) k B T/R co . 

To illustrate the difference between a and Xf, we in- 
clude in Fig. [12] the two asymptotical power spectra aris- 
ing from Af as the dash-dotted line, and from a when 
simply inserted into Eq. (|4Tj) as the dashed line. Xf nicely 
describes the measured fluctuation spectrum, whereas a 
results in a significantly damped spectrum. Thus, there is 
a notable difference between a and A^ of approximately 
30 % even for the small system size considered in the 
simulation. 

Finally, we note that the line tensions, A, of membranes 
formed by biologically relevant lipids are typically on the 
order of 10 pN i 19 ' 22 ' 100 Using the previously determined 
scale factors for length and energy, we find that 1 pN ~ 
0.5 kBT/R co , i.e. both simulation results A « 10(1) pN 
and Xf ~ 11(1) pN match the experimental values. 

VI. CONCLUSIONS 

In this study, we have presented a soft, solvent-free 
coarse-grained model for the simulation of lipid bilay- 
ers. The non-bonded interactions are inspired by a field- 
theoretic description and take the form of a third-order 
expansion of the excess free energy functional in the den- 
sities of the hydrophilic and hydrophobic beads. The 
numerical values of the expansion coefficients are related 
to a few thermodynamic key characteristics like the den- 
sity of the hydrophobic core and its compressibility. 30 
The local structural properties of the bilayer, i.e., liquid- 
like packing of the coarse-grained beads, can be indepen- 
dently adjusted by means of weighting functions for the 
attractive pair-wise interactions and the repulsive triple 
interactions. In this way, we devise a flexible model 
where the interaction parameters bear a clear physical 



interpretation. With the described DPD simulation tech- 
nique, it is possible to simulate bilayer patches of a size 
of A = 100 x 100 dih 2 for up to At = 1 ms. 

Depending on the molecular asymmetry, the lipids self- 
assemble into spherical and cylindrical micelles, wormlike 
micelles, bilayers, and inverted structures. Static and 
dynamic properties of the bilayer membranes have been 
calculated. In the fluid phase, we obtained a bending 
rigidity of k ~ 15 — 19 fcsT, an area compression modulus 
of Ha ~ 70 mN/m, and a molecular aspect ratio of 1 : 7. 
These numbers match the orders of magnitude observed 
in experiments. 

In the second part, the phase behavior of the bilayer 
has been investigated. Depending on the harshness of the 
short-range repulsive interactions and on the density of 
the hydrophobic interior, we observe the fluid phase, L a , 
and three different gel phases. We have studied the main 
phase transition, Lp ■<-> L a , in detail. By means of Um- 
brella Sampling and the Weighted Histogram Analysis 
Method, we have calculated the free energy as a function 
of a conformational order parameter near the phase co- 
existence and obtained the free energy profile across the 
miscibility gap. 

The phase coexistence has been accurately located by 
three different, although not independent, methods giv- 
ing consistent results: (i) histogram reweighting, (ii) free 
energy perturbation calculations, and (iii) the calculation 
of the specific heat. Different geometries of the minor- 
ity phase, like droplets and slabs, have been observed in 
the miscibility gap. Finally, the line tension separating 
different domains has been calculated by means of two 
different methods, giving a final result of A m 10 pN, 
which is in excellent agreement with experimental stud- 
ies. The computational method outline in the present 
work is not restricted to studying phase transitions of 
single-component lipid membranes. For instance, one 
could study fluid-fluid coexistence in lipid mixtures or 
asymmetric bilayers. In such a case, one will setup simu- 
lations at constant lateral tension in a semigrand canon- 
ical ensemble, where a suitable order parameter will be 
the overall lipid composition. 

We hope that the soft, coarse-grained model for lipid 
bilayer membranes and the computational techniques will 
find further applications in the study of collective phe- 
nomena in membranes. 
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Appendix A: Integration algorithm 



with P t = {<jyy + a zz )/2, where 



AL X / — 1 \mA J 



Most of the simulations have been performed in a sta- 
tistical ensemble, in which the average tangential pres- 
sure Pt and the height L x of the simulation box are 
kept constant, so that the area A(t) enters as one addi- 
tional dynamic degree of freedom. Although several re- 
lated algorithms have been published ) 101 ' 102 we describe 
our symplectic integration algorithms for completeness 
in this appendix. It is derived from the Langevin 
piston method developed in Ref. 64 A second integra- 
tion algorithm with two additional degrees of freedom 
(L y (t), L z (t)) can be derived in a similar way. It is used 
to simulate bilayers in the gel phase, where isotropic fluc- 
tuations of the area are inappropriate due to the hexag- 
onal ordering of the lipids. For brevity the details of this 
similar algorithm are omitted. 

In a thermally isolated system the first law of thermo- 
dynamics reads 



dE = -PdV + -/dA = 
= -F\L x dA, 



-PdV + L X (P - P t )dA 



where P is the normal pressure and 7 = L X (P — P t ) is 
the surface tension. Hence, dH = d(E + P t L x A) = 
and the enthalpy H = E + PtL x A is a conserved quan- 
tity. Following Anderson^22 we now introduce scaled co- 
ordinates jji — Si y \fA and Zi — Si Z \fA tangential to the 
plane, retaining the normal coordinates, Xi. In the parti- 
cle velocities, iji = hi y \J~A the second term is deliberately 
omitted to achieve independent fluctuations of A and Si y . 
One can now write down the classical Lagrangian, 



2 + Asl 



Asl) 



QLl 



A 



-U(A,{s }) - P t L x A 



where we have introduced an artificial mass Q for the 
new degree of freedom, setting the timescale of the area's 
fluctuations. Introducing the canonically conjugated mo- 
menta pi X — mx,Hi a — mAsi a ,TTA — QL X A, with 
a = {y,z}, one can derive the Hamiltonian by means 
of a Legendre transformation 



/ Pi 



^ \ 2m 2mA 2mA j 

+U (A,{s iy ,s iz })+T t L x A. 
The Hamilton equations of motion read 

F 



2QLI 



Pix 

m 



Pi. 



= TT. la = F ia VA 

mA 



A 



TT A = L x (P t - 7\) 



are the diagonal entries in the pressure tensor, that can 
be calculated with the virial theorem. 

Following Tuckerman et al.^2^, we split the Liouville 
operator, iL, into a sum of simpler operators 

Odd 

Fix o + Fiy VA- — + F iz VA 



d 
dir A 



dp lx d-Ki. 

iL 2 = L x (P t -P7) 

f - nA _ 
3 " QLl OA 

.f Tr^Pix d t ir iy d : ir. lz d 

iL 4 = 2^ 



dm 



m dxi mA dsi y mA d-Si 



We approximate the unitary time evolution operator by 
the Trotter factorization, yielding 



U{At) 



_ g iLAt _ e i(L 1 +L 2 +L 3 +L 4 )At 



Applying these operators one after another from the right 
to the left onto the phase space vector, one obtains a 
symplectic integration algorithm in the canonically con- 
jugated quantities: 



7Tia(-jr) 


= p ix (0) + 

= 7T m (0)H 


F- At. 

1 ix 2 

-F iay /A(0)f 




= 7T A (0) + 


L x (P t -p- t )f 


A(f) = 


- A(0) + ^<£'/ 2) A 2 * 


Xi(At) = 
s fa (At) 


= x l (0) + £ 
= s ia (0) + 


m 
mA 


A(At) = 


--A(At) + 


TT A (At/2) At 
QLl 2 


7T A {At) 


= Mf)- 


+• L x {P t -F\)f 


Pix(At) 

TTi a (At) 


= Pi*(f) 

-7T- (At) 


+ F -At 

\ J - %x 2 

+ F ia ^A(At)f 



To simplify the usage of this algorithm the scaled coordi- 
nates, Si a , are finally substituted by the real coordinates, 
i.e. 

Vi (t)=s iv (t)^A(f} Piy{t) = ^K 



Zi(t) = »i,(t)VAW Pi*(t) = 



This invokes some rescaling steps in the final algorithm. 
Hence: 
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1. Calculation of the temporary momenta using the 
old forces: 



Pic 

(At/2) = Plx (0) + F l 



At 



= n la (At/2) _ At 



2. Calculation of the intermediate area momentum us- 
ing the tangential pressure Pt evaluated with the 
old forces F(0) and the new temporary momenta 



Pi 



n A (At/2) = tt a (0) + L x (P t - P t ) 



At 



3. First half of the integration of the area: 

4/ « i , -K A (Atl2) At 

A(At/2) = A(0) + M Q J ' — 

4. Integration of the particle coordinates: 

Xi(At) = x l (0) + Ptx{At/2) At 
m 

a> , s ia (At)VW) = a i (0) + ^ ) P ^t 

5. Second half of the integration of the area: 

QLi 2 

6. Rescaling of the coordinates according to: 



7. Recalculation of the new forces using the new co- 
ordinates and recalculation of the pressure tensor 
using the new forces and the temporary momenta 
p" ■ 

8. Calculation of the final area momentum using the 
tangential pressure Pt evaluated with the new 
forces F(At) and the temporary momenta p'/ a : 



ir A (At) = tta(A*/2) + L x (P t - P t ) 



At 



9. Final integration of the particle's momentum: 
Pix {At) = p tx {At/2) + F ix ^- 

T^ia (At) = P" a + F ia ^ 



Up to now the integration algorithm has been formu- 
lated in the microcanonical ensemble, where the total 
energy is conserved. The switch to the iV Pi T-ensemble 
is performed with the DPD thermostat for the particle 
interactions and with a Langevin thermostat for the area 
A(t). The former is described in section HTl and the latter 
is achieved via the replacement^ 



L x {P t 



— x At 



L x (P t 



-1a 



At 



QLI 2 



+ ^k B T lA At/2U- 



In the last equation, £ A is a random number drawn from 
a uniform distribution with (£4) = and (£ A ) = 1 and 
j A is a friction coefficient. 

In our simulations we have used the values Q = 
0.0001,7a = 0.1 and Pt = 0, which corresponds to a 
simulation at vanishing lateral tension. 



Appendix B: Capillary waves and line tension 

In this Appendix we briefly present the derivation of 
Eqs. ([41~j) and (|43|. The statistical properties of the phase 
boundary follow from the capillary wave Hamiltonian 



X J d.r ] 





[h'(x)y 



XL+^jdx [h'(x)}\ 




where A is the bare line tension, L denotes the projected 
length of the line, and h(x) represents the position of 
the phase boundary. The coordinate system is chosen 
such that the mean position of the boundary vanishes, 
i.e. (h(x)) — 0. Expanding 



h(x) = with 1 



2irn 



in a Fourier series with wave numbers q, n € Z, that are 
commensurate with the periodic boundary conditions, 
%cap becomes diagonal in g-space and the different modes 
decouple. Hence, we rewrite "H C ap in the form 



\L + Uf with H f 



XL 



2 2 

91 q ■ 



The Hamiltonian fif is the starting point for all further 
calculations. On the one hand, the fluctuation power 
spectrum is readily obtained from this expression using 
the equipartition theorem^ 



XLq 2 



On the other hand, we calculate the free energy con- 
tribution, AF = F — XL, from the fluctuations of the 
boundary. The canonical partition function, Z, involves 



19 



a functional integral over all possible interface profiles, 
h(x), which is equivalent to a functional integral over all 
complex Fourier coefficients, h q : 



J V [h q ] cxp (- 



n ^iexp 



XL 
2k B T 



Since h(x) is a real-valued function, the complex coeffi- 
cients h q possess the Hermitian redundancy, i.e. h— q = 
h q , so that we decompose the functional integral into sep- 
arate integrals over the real and the imaginary parts of 
h q using only the positive modes, q > 0. Thus, 



n/%=n 



g>0 



d (Re h q ) d (lmh q ) 



We evaluate Z by carrying out the Gaussian quadratures, 
yielding 



2irk R T 1 



q>0 



XL 



To compute the free energy AF = —k B T\nZ, the range 
of wave vectors has to be restricted. The simplest method 
is to introduce cut-offs for small q as well as for high 
q. The small q-cutoff naturally arises from the periodic 



2tt/L. The UV-cutoff 
= L/a is more difficult 



boundary conditions, i.e., q m \ n = 

?max 27Tn nlSL ^f L 1lT j CI, Tl max 

because it introduces an additional length scale. This 
length scale characterizes the smallest length scale, on 
which the fluctuations of the boundary line can be de- 
scribed by the capillary wave Hamiltonian. Using these 
two cut-offs, we obtain 



= - hxZ = 2^ In ( - , - • 4ttV 



k B T 



In 



n=l 

2ttXL 
k R T 



2nk B T 
f 21n(ra max !). 



AF is an extensive quantity, proportional to the size, L, 
of the system: 

AF _ _ lf (2ttXL 3 \ , C 

LkB~r- a [ \yk^f)~ J = ^T 

Hence, the total interfacial free energy is given by 

F = XL + AF = (A + C)L = oL. 

Thus, the bare line tension, A, in a one-dimensional sys- 
tem differs from the thermodynamic line tension, a, by 
a quantity, C, which stems from the fluctuations of the 
contact line and which logarithmically depends on L and 



* Electronic address: mmueller@theorie.physik.uni-goettingen.de 14 

1 H. Lodish, A. Berk, C. A. Kaiser, M. Krieger, M. P. 
Scott, and A. Bretscher, Molecular Cell Biology (Palgrave 1 
Macmillan, New York, 2007). 16 

2 O. G. Mouritsen, Life - As a Matter of Fat (Springer, 
Berlin, 2005). 17 

3 M. Miiller, K. Katsov, and M. Schick, J. Polym. Sci. B: 
Polymer Physics 41, 1441 (2003). 

4 M. Muller, K. Katsov, and M. Schick, Phys. Rep. 434, 18 
113 (2006). 

5 G. Brannigan, L. C. L. Lin, and F. L. H. Brown, Eur. 19 
Biophys. J. 35, 104 (2006). 

6 M. Venturoli, M. M. Sperotto, M. Kranenburg, and 20 
B. Smit, Phys. Rep. 437, 1 (2006). 

7 R. Koynova and M. Caffrey, Biochim. Biophys. Acta Rev. 
Biomembranes 1376, 91 (1998). 

8 J. F. Nagle and S. Tristram-Nagle, Biochim. Biophys. 22 
Acta Rev. Biomembranes 1496, 159 (2000). 

9 J. H. Ipsen, K. Jorgensen, and O. G. Mouritsen, Biophys. 
J. 58, 1099 (1990). 

10 S. J. Marrink and A. E. Mark, Biophys. J. 87, 3894 24 

(2004) . 

11 S. J. Marrink, J. Risselada, and A. E. Mark, Chem. Phys. 25 
Lip. 135, 223 (2005). 

12 M. Kranenburg and B. Smit, J. Phys. Chem. B 109, 6553 26 

(2005) . 

13 O. Lenz and F. Schmid, J. Mol. Liquids 117, 147 (2005). 27 



I. R. Cooke and M. Deserno, J. Chem. Phys. 123, 224710 

(2005) . 

M. J. Stevens, J. Chem. Phys. 121, 11942 (2004). 

J. D. Revalee, M. Laradji, and P. B. Sunil Kumar, J. 

Chem. Phys. 128, 035102 (2008). 

O. G. Mouritsen, A. Boothroyd, R. Harris, N. Jan, 

T. Lookman, L. MacDonald, D. A. Pink, and M. J. Zuck- 

ermann, J. Chem. Phys. 79, 2027 (1983). 

J.-M. Allain, C. Storm, A. Roux, M. B. Amar, and J.-F. 

Joanny, Phys. Rev. Lett. 93, 158104 (2004). 

Joannis, F. Y. Jiang, and J. T. Kindt, Langmuir 22, 998 

(2006) . 

T. Baumgart, S. T. Hess, and W. W. Webb, Nature 425, 
821 (2003). 

A. Tian, C. Johnson, W. Wang, and T. Baumgart, Phys. 
Rev. Lett. 98, 208102 (2007). 

C. Esposito, A. Tian, S. Melamed, C. Johnson, S.-Y. Tee, 

and T. Baumgart, Biophys. J. 93, 3169 (2007). 

G. M. Torrie and J. P. Valleau, J. Comput. Phys. 23, 187 

(1977). 

A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 
63, 1195 (1989). 

P. Virnau and M. Muller, J. Chem. Phys. 120, 10925 
(2004). 

A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 
61, 2635 (1988). 

S. Kumar, D. Bouzida, R. Swendsen, P. Kollman, and 



20 



J. Rosenberg, J. Comp. Chem. 13, 1011 (1992). 

28 M. Souaille and B. Roux, Comp. Phys. Comm. 135, 40 
(2001). 

29 C. Chipot and A. Pohorille, eds., Free energy calcula- 
tions: theory and applications in chemistry and biology 
(Springer, Berlin, 2007). 

30 K. Ch. Daoulas and M. Miiller, Adv. Polym. Sci 224, 197 
(2009). 

31 K. Ch. Daoulas and M. Miiller, J. Chem. Phys. 125, 
184904 (2006). 

32 L. Gao, J. Shillcock, and R. Lipowsky, J. Chem. Phys. 
126, 015101 (2007). 

33 J. M. Drouffe, A. C. Maggs, and S. Leibler, Science 254, 
1353 (1991). 

34 H. Noguchi and M. Takasu, Phys. Rev. E 64, 041913 
(2001). 

35 O. Farago, J. Chem. Phys. 119, 596 (2003). 

36 G. Brannigan and F. L. H. Brown, J. Chem. Phys. 120, 
1059 (2004). 

37 A. A. Louis, J. Phys.: Condens. Matter 14, 9187 (2002). 

38 M. Miiller and L. G. MacDowell, Macromolecules 33, 3902 

(2000) . 

39 M. Miiller, L. G. MacDowell, P. Virnau, and K. Binder, 
J. Chem. Phys. 117, 5480 (2002). 

40 P. H. van Konynenburg and R. L. Scott, Phil. Trans. R. 
Soc. A 298, 496 (1980). 

41 J. Wang and M. Miiller, Macromolecules 42, 2251 (2009). 

42 J. F. Wang and M. Miiller, J. Phys. Chem. B 113, 11384 
(2009). 

43 T. S. Khasanshin, A. P. Shchamialiou, and O. G. Pod- 
dubskij, Int. J. Thermophys. 24, 1277 (2003). 

44 O. G. Mouritsen, Chem. Phys. Lip. 57, 179 (1991). 

45 W. A. Curtin and N. W. Ashcroft, Phys. Rev. A 32, 2909 
(1985). 

46 F. van Swol and J. R. Henderson, Phys. Rev. A 43, 2932 
(1991). 

47 A. Yethiraj, J. Chem. Phys. 109, 3269 (1998). 

48 M. Miiller, L. G. MacDowell, and A. Yethiraj, J. Chem. 
Phys. 118, 2929 (2003). 

49 J. D. Weeks, K. Katsov, and K. Vollmayr, Phys. Rev. 
Lett. 81, 4400 (1998). 

50 K. Katsov and J. D. Weeks, J. Phys. Chem. B 105, 6738 

(2001) . 

51 P. Warren and P. Espanol, Europhys. Lett 30, 191196 
(1995). 

52 C. N. Likos, B. M. Mladek, D. Gottwald, and G. Kahl, J. 
Chem. Phys. 126, 224502 (2007). 

53 B. M. Mladek, M. J. Fernaud, G. Kahl, and M. Neumann, 
Condensed Matter Physics 8, 135 (2005). 

54 B. M. Mladek, D. Gottwald, G. Kahl, M. Neumann, and 
C. N. Likos, J. Phys. Chem. B 111, 12799 (2007). 

55 W. A. Curtin and N. W. Ashcroft, Phys. Rev. A 32, 2909 
(1985). 

56 I. Pagonabarraga and D. Frenkel, J. Chem. Phys. 115, 
5015 (2001). 

57 S. Y. Trofimov, E. L. F. Nies, and M. A. J. Michels, J. 
Chem. Phys. 117, 9383 (2002). 

58 P. B. Warren, Phys. Rev. E 68, 066702 (2003). 

59 P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. 
Lett 19, 155 (1992). 

60 J. M. V. A. Koelman and P. J. Hoogerbrugge, Europhys. 
Lett 21, 363 (1993). 

61 B. Diinweg and W. Paul, Int. J. Mod. Phys. C 2, 817 
(1991). 



62 P. Espanol and P. Warren, Europhys. Lett. 30, 191 (1995). 

63 M. P. Allen and D. J. Tildesley, Computer Simulation of 
Liquids (Clarendon Press, Oxford, 1987). 

64 A. Kolb and B. Diinweg, J. Chem. Phys. Ill, 4453 (1999). 

65 S. Plimpton, J. Comp. Phys. 117, 1 (1995). 

66 J. N. Israelachvili, Intermolecular and Surfaces Forces 
(Academic Press, London, 1991), 2nd ed. 

67 J. F. Nagle, R. Zhang, S. Tristram-Nagle, W. Sun, H. I. 
Petrache, and R. M. Suter, Biophys. J. 70, 1419 (1996). 

68 H. I. Petrache, S. W. Dodd, and M. F. Brown, Biophys. 
J. 79, 3172 (2000). 

69 W. K. den Otter, J. Chem. Phys. 123, 214906 (2005). 

70 An example of such a configuration is depicted in Fig. [9)3. 

71 E. Lindahl and O. Edholm, Biophys. J. 76, 426 (2000). 

72 D. Marsh, Chem. Phys. Lipids 144, 146 (2006). 

73 P. B. Canham, J. Theor. Bio. 26, 61 (1970). 

74 W. Helfrich, Zeitschrift Naturforschung C 28, 693 (1973). 

75 E. A. Evans, Biophys. J. 14, 923 (1974). 

76 M. Miiller and M. Schick, J. Chem. Phys. 105, 8885 
(1996). 

77 E. Boek, J. Padding, W. den Otter, and W. Briels, J. 
Phys. Chem. B 109, 19851 (2005). 

78 G. Brannigan and F. L. H. Brown, Biophys. J. 90, 1501 

(2006) . 

79 U. Seifert, Adv. Phys. 46, 13 (1997). 

80 R. Goetz, G. Gompper, and R. Lipowsky, Phys. Rev. Lett. 
82, 221 (1999). 

81 E. Reister and U. Seifert, Europhys. Lett 71, 859 (2005). 

82 E. Reister-Gottfried, S. M. Leitenberger, and U. Seifert, 
Phys. Rev. E 75, 011908 (2007). 

83 N. Kahya, D. Scherfeld, K. Bacia, and P. Schwille, J. 
Struct. Biol. 147, 77 (2004). 

84 M. S. C. Abreu, M. Joao Moreno, and W. L. C. Vaz, 
Biophys. J. 87, 353 (2004). 

85 J. Liu, S. Qi, J. Groves, and A. Chakraborty, J. Phys. 
Chem. B 109, 199960 (2005). 

86 S. Leekumjorn and A. K. Sum, Biochim. Biophys. Acta 
Biomembranes 1768, 354 (2006). 

87 K. Binder, Phys. Rev. A 25, 1699 (1982). 

88 B. A. Berg, U. Hansmann, and T. Neuhaus, Z. Phys. B 
90, 229 (1993). 

89 J. E. Hunter and W. P. Reinhardt, J. Chem. Phys. 103, 
8627 (1995). 

90 K. Binder, Physica A-Statistical Mechanics and Its Ap- 
plications 319, 99 (2003). 

91 L. G. MacDowell, P. Virnau, M. Miiller, and K. Binder, 
J. Chem. Phys. 120, 5293 (2004). 

92 D. S. Dean and R. R. Horgan, Phys. Rev. E 76, 041102 

(2007) . 

93 R. Golestanian, Europhys. Lett 36, 557 (1996). 

94 M. P. Gelfand and M. E. Fisher, Physica. A 166, 1 (1990). 

95 S. A. Safran, Statistical thermodynamics of surfaces, in- 
terfaces and membranes (Addison Wesley, Reading MA, 
1994). 

96 E. Chacon and P. Tarazona, J. Phys.: Condens. Matter 
17, S3493 (2005). 

97 C. Pastorino, K. Binder, and M. Miiller, Macromolecules 
42, 401 (2009). 

98 A. Werner, F. Schmid, M. Miiller, and K. Binder, Phys. 
Rev. E 59, 728 (1999). 

99 M. Miiller and J. J. de Pablo, Lec. Notes Phys. 703, 67 
(2006). 

100 E. Karatekin, O. Sandre, H. Guitouni, N. Borghi, P. H. 
Puech, and F. Brochard-Wyart, Biophys. J. 84, 1734 



21 



(2003). 

A. F. Jakobsen, O. G. Mouritsen, and G. Besold, J. Chem. 
Phys. 122, 204901 (2005). 

Y. Zhang, S. E. Feller, B. R. Brooks, and R. W. Pastor, 

J. Chem. Phys. 103, 10252 (1995). 

H. Anderson, J. Chem. Phys. 72, 2384 (1980). 

M. Tuckcrman, B. Berne, and G. Martyna, J. Chem. 



Phys. 97, 1990 (1992). 

15 W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph- 
ics 14, 33 (1996). 

16 W. Shinoda and S. Okazaki, J. Chem. Phys. 109, 1517 
(1998). 



